ATM 623: Climate Modeling¶
Brian E. J. Rose, University at Albany
Lecture 7: Grey radiation modeling with climlab¶
About these notes:¶
This document uses the interactive IPython notebook format (now also called Jupyter). The notes can be accessed in several different ways:
- The interactive notebooks are hosted on
githubat https://github.com/brian-rose/ClimateModeling_courseware - The latest versions can be viewed as static web pages rendered on nbviewer
- A complete snapshot of the notes as of May 2015 (end of spring semester) are available on Brian's website.
Many of these notes make use of the climlab package, available at https://github.com/brian-rose/climlab
Contents¶
- Introducing
climlab - Using
climlabto implement the two-layer leaky greenhouse model - The observed annual, global mean temperature profile
- A 30-layer model using the observed temperatures
- Radiative forcing in the 30-layer model
- Radiative equilibrium in the 30-layer model
- Radiative-Convective Equilibrium in the 30-layer model
- Putting stratospheric ozone in the grey-gas model
1. Introducing climlab¶
climlab is a flexible engine for process-oriented climate modeling.
It is based on a very general concept of a model as a collection of individual,
interacting processes. climlab defines a base class called Process, which
can contain an arbitrarily complex tree of sub-processes (each also some
sub-class of Process). Every climate process (radiative, dynamical,
physical, turbulent, convective, chemical, etc.) can be simulated as a stand-alone
process model given appropriate input, or as a sub-process of a more complex model.
New classes of model can easily be defined and run interactively by putting together an
appropriate collection of sub-processes.
climlab is a work-in-progress, and the code base will evolve substantially over the course of this semester.
The latest code can always be found on github:
https://github.com/brian-rose/climlab
You are strongly encouraged to clone the climlab repository and use git to keep your local copy up-to-date.
Running this notebook requires that climlab is already installed on your system.
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import netCDF4 as nc
import climlab
2. Using climlab to implement the two-layer leaky greenhouse model¶
One of the things that climlab is set up to do is the grey-radiation modeling we have already been discussing.
Since we already derived a complete analytical solution to the two-layer leaky greenhouse model, we will use this to validate the climlab code.
Validation¶
We want to verify that the model reproduces the observed OLR given observed temperatures, and the absorptivity that we tuned in the analytical model. The target numbers are:
\begin{align} T_s &= 288 \text{ K} \\ T_0 &= 275 \text{ K} \\ T_1 &= 230 \text{ K} \\ \end{align}$$ \epsilon = 0.58377 $$$$ OLR = 239 \text{ W m}^{-2} $$Initialize a model in climlab¶
The first thing we do is create a new model.
The following example code is sparsely commented but will hopefully orient you on the basics of defining and working with a climlab Process object.
# Test in a 2-layer atmosphere
col = climlab.GreyRadiationModel(num_lev=2)
print col
col.subprocess
col.state
col.Ts
col.Ts[:] = 288.
col.Tatm[:] = np.array([275., 230.])
col.state
LW = col.subprocess['LW']
print LW
LW.absorptivity
LW.absorptivity = 0.58377
LW.absorptivity
# This does all the calculations that would be performed at each time step,
# but doesn't actually update the temperatures
col.compute_diagnostics()
col.diagnostics
# Check OLR against our analytical solution
col.diagnostics['OLR']
col.state
# perform a single time step
col.step_forward()
col.state
# integrate out to radiative equilibrium
col.integrate_years(2.)
# Check for equilibrium
col.diagnostics['ASR'] - col.diagnostics['OLR']
# Compare these temperatures against our analytical solutions for radiative equilibrium
col.state
So it looks like climlab agrees with our analytical results. That's good.
3. The observed annual, global mean temperature profile¶
We want to model the OLR in a column whose temperatures match observations. As we've done before, we'll calculate the global, annual mean air temperature from the NCEP Reanalysis data.
# This will try to read the data over the internet.
# If you have a local copy of the data, just use the local path to the .nc file instead of the URL
ncep_url = "http://www.esrl.noaa.gov/psd/thredds/dodsC/Datasets/ncep.reanalysis.derived/"
ncep_air = nc.Dataset( ncep_url + "pressure/air.mon.1981-2010.ltm.nc" )
level = ncep_air.variables['level'][:]
lat = ncep_air.variables['lat'][:]
# A log-pressure height coordinate
zstar = -np.log(level/1000)
Tzon = np.mean(ncep_air.variables['air'][:],axis=(0,3))
Tglobal = np.average( Tzon , weights=np.cos(np.deg2rad(lat)), axis=1) + climlab.constants.tempCtoK
# Note the useful conversion factor. climlab.constants has lots of commonly used constant pre-defined
# Here we are plotting with respect to log(pressure) but labeling the axis in pressure units
fig = plt.figure( figsize=(8,6) )
ax = fig.add_subplot(111)
ax.plot( Tglobal , zstar )
yticks = np.array([1000., 750., 500., 250., 100., 50., 20., 10.])
ax.set_yticks(-np.log(yticks/1000.))
ax.set_yticklabels(yticks)
ax.set_xlabel('Temperature (K)', fontsize=16)
ax.set_ylabel('Pressure (hPa)', fontsize=16 )
ax.set_title('Global, annual mean sounding from NCEP Reanalysis', fontsize = 24)
ax.grid()
# initialize a grey radiation model with 30 levels
col = climlab.GreyRadiationModel()
print col
# interpolate to 30 evenly spaced pressure levels
lev = col.lev
Tinterp = np.flipud(np.interp(np.flipud(lev), np.flipud(level), np.flipud(Tglobal)))
Tinterp
# Need to 'flipud' because the interpolation routine needs the pressure data to be in increasing order
# Initialize model with observed temperatures
col.Ts[:] = Tglobal[0]
col.Tatm[:] = Tinterp
# A handy re-usable routine for making a plot of the temperature profiles
# We will plot temperatures with respect to log(pressure) to get a height-like coordinate
def plot_sounding(collist):
color_cycle=['r', 'g', 'b', 'y']
# col is either a column model object or a list of column model objects
if isinstance(collist, climlab.Process):
# make a list with a single item
collist = [collist]
fig = plt.figure()
ax = fig.add_subplot(111)
for i, col in enumerate(collist):
zstar = -np.log(col.lev/climlab.constants.ps)
ax.plot(col.Tatm, zstar, color=color_cycle[i])
ax.plot(col.Ts, 0, 'o', markersize=12, color=color_cycle[i])
#ax.invert_yaxis()
yticks = np.array([1000., 750., 500., 250., 100., 50., 20., 10.])
ax.set_yticks(-np.log(yticks/1000.))
ax.set_yticklabels(yticks)
ax.set_xlabel('Temperature (K)')
ax.set_ylabel('Pressure (hPa)')
ax.grid()
return ax
# This should look just like the observations
plot_sounding(col)
Tune absorptivity to get observed OLR¶
col.compute_diagnostics()
col.diagnostics['OLR']
# Need to tune absorptivity to get OLR = 239
epsarray = np.linspace(0.01, 0.1, 100)
OLRarray = np.zeros_like(epsarray)
for i in range(epsarray.size):
col.subprocess['LW'].absorptivity = epsarray[i]
col.compute_diagnostics()
OLRarray[i] = col.diagnostics['OLR']
plt.plot(epsarray, OLRarray)
plt.grid()
The necessary value seems to lie near 0.055 or so.
We can be more precise with a numerical root-finder.
def OLRanom(eps):
col.subprocess['LW'].absorptivity = eps
col.compute_diagnostics()
return col.diagnostics['OLR'] - 239.
# Use numerical root-finding to get the equilibria
from scipy.optimize import brentq
# brentq is a root-finding function
# Need to give it a function and two end-points
# It will look for a zero of the function between those end-points
eps = brentq(OLRanom, 0.01, 0.1)
print eps
col.subprocess['LW'].absorptivity = eps
col.subprocess['LW'].absorptivity
col.compute_diagnostics()
col.diagnostics['OLR']
5. Radiative forcing in the 30-layer model¶
Let's compute radiative forcing for a 2% increase in absorptivity.
# clone our model using a built-in climlab function
col2 = climlab.process_like(col)
print col2
col2.subprocess['LW'].absorptivity *= 1.02
col2.subprocess['LW'].absorptivity
# Radiative forcing by definition is the change in TOA radiative flux, HOLDING THE TEMPERATURES FIXED.
col2.Ts - col.Ts
col2.Tatm - col.Tatm
col2.compute_diagnostics()
col2.diagnostics['OLR']
The OLR decreased after we added the extra absorbers, as we expect. Now we can calculate the Radiative Forcing:
RF = -(col2.diagnostics['OLR'] - col.diagnostics['OLR'])
print 'The radiative forcing is %f W/m2.' %RF
re = climlab.process_like(col)
# To get to equilibrium, we just time-step the model forward long enough
re.integrate_years(2.)
# Check for energy balance
re.diagnostics['ASR'] - re.diagnostics['OLR']
plot_sounding([col, re])
Some properties of the radiative equilibrium temperature profile:
- The surface is warmer than observed.
- The lower troposphere is colder than observed.
- Very cold air is sitting immediately above the warm surface.
- There is no tropopause, no stratosphere.
7. Radiative-Convective Equilibrium in the 30-layer model¶
We recognize that the large drop in temperature just above the surface is unphysical. Parcels of air in direct contact with the ground will be warmed by mechansisms other than radiative transfer.
These warm air parcels will then become buoyant, and will convect upward, mixing their heat content with the environment.
We parameterize the statistical effects of this mixing through a convective adjustment.
At each timestep, our model checks for any locations at which the lapse rate exceeds some threshold. Unstable layers are removed through an energy-conserving mixing formula.
This process is assumed to be fast relative to radiative heating. In the model, it is instantaneous.
rce = climlab.RadiativeConvectiveModel(adj_lapse_rate=6.)
print rce
This model is exactly like our previous models, except for one additional subprocess called convective adjustment.
We passed a parameter adj_lapse_rate (in K / km) that sets the neutrally stable lapse rate -- in this case, 6 K / km.
This number is chosed to very loosely represent the net effect of moist convection. We'll look at this in more detail later.
# Set our tuned absorptivity value
rce.subprocess['LW'].absorptivity = eps
# Run out to equilibrium
rce.integrate_years(2.)
# Check for energy balance
rce.diagnostics['ASR'] - rce.diagnostics['OLR']
# Make a plot to compare observations, Radiative Equilibrium, and Radiative-Convective Equilibrium
plot_sounding([col, re, rce])
Introducing convective adjustment into the model cools the surface quite a bit (compared to Radiative Equilibrium, in green here) -- and warms the lower troposphere. It gives us a MUCH better fit to observations.
But of course we still have no stratosphere.
8. Putting stratospheric ozone in the grey-gas model¶
Our model has no equivalent of the stratosphere, where temperature increases with height. That's because our model has been completely transparent to shortwave radiation up until now.
We can load the observed ozone climatology from the input files for the CESM model:
datapath = "http://ramadda.atmos.albany.edu:8080/repository/opendap/latest/Top/Users/Brian+Rose/CESM+runs/"
endstr = "/entry.das"
ozone = nc.Dataset( datapath + 'som_input/ozone_1.9x2.5_L26_2000clim_c091112.nc' + endstr )
print ozone.variables['O3']
lat_O3 = ozone.variables['lat'][:]
lon_O3 = ozone.variables['lon'][:]
lev_O3 = ozone.variables['lev'][:]
The pressure levels in this dataset are:
print lev_O3
Take the global average of the ozone climatology, and plot it as a function of pressure (or height)¶
O3_zon = np.mean( ozone.variables['O3'][:],axis=(0,3) )
print O3_zon.shape
O3_global = np.average( O3_zon, axis=1, weights=np.cos(np.deg2rad(lat_O3)))
print O3_global
ax = plt.figure(figsize=(10,8)).add_subplot(111)
ax.plot( O3_global * 1.E6, -np.log(lev_O3/climlab.constants.ps) )
ax.set_xlabel('Ozone (ppm)', fontsize=16)
ax.set_ylabel('Pressure (hPa)', fontsize=16 )
ax.set_yticks( -np.log(yticks/1000.) )
ax.set_yticklabels( yticks )
ax.grid()
ax.set_title('Global, annual mean ozone concentration', fontsize = 24);
This shows that most of the ozone is indeed in the stratosphere, and peaks near the top of the stratosphere.
Now create a new column model object on the same pressure levels as the ozone data. We are also going set an adjusted lapse rate of 6 K / km.
oz_col = climlab.RadiativeConvectiveModel(lev = lev_O3, adj_lapse_rate=6)
print oz_col
Now we will do something new: let the column absorb some shortwave radiation. We will assume that the shortwave absorptivity is proportional to the ozone concentration we plotted above.
First we have to deal with a little inconsistency:
print lev_O3
print oz_col.lev
The two arrays are in reverse order!
So we need to flip the ozone data before using it:
O3_flipped = np.flipud(O3_global)
Now we need to weight the absorptivity by the pressure (mass) of each layer.
# This number is an arbitrary parameter that scales how absorptive we are making the ozone
# in our grey gas model
ozonefactor = 75
dp = oz_col.Tatm.domain.lev.delta
epsSW = np.flipud(O3_global) * dp * ozonefactor
We want to use the field epsSW as the absorptivity for our SW radiation model.
Let's see what the absorptivity is current set to:
print oz_col.subprocess['SW'].absorptivity
It defaults to zero.
Before changing this (putting in the ozone), let's take a look at the shortwave absorption in the column:
oz_col.compute_diagnostics()
oz_col.diagnostics['SW_absorbed_atm']
Let's now put in the ozone:
oz_col.subprocess['SW'].absorptivity = epsSW
print oz_col.subprocess['SW'].absorptivity
Let's check how this changes the SW absorption:
oz_col.compute_diagnostics()
oz_col.diagnostics['SW_absorbed_atm']
It is now non-zero, and largest near the top of the column (bottom of array) where the ozone concentration is highest.
Now it's time to run the model out to radiative-convective equilibrium
oz_col.integrate_years(1.)
print oz_col.diagnostics['ASR'] - oz_col.diagnostics['OLR']
And let's now see what we got!
# Make a plot to compare observations, Radiative Equilibrium, Radiative-Convective Equilibrium, and RCE with ozone!
plot_sounding([col, re, rce, oz_col])
And we finally have something that looks looks like the tropopause, with temperature increasing above at approximately the correct rate.
There are still plenty of discrepancies between this model solution and the observations, including:
- Tropopause temperature is too warm, by about 15 degrees.
- Surface temperature is too cold
There are a number of parameters we might adjust if we wanted to improve the fit, including:
- Longwave absorptivity
- Surface albedo
Feel free to experiment! (That's what models are for, after all).
The take home message¶
The dominant effect of stratospheric ozone is to vastly increase the radiative equilibrium temperature in the ozone layer. The temperature needs to be higher so that the longwave emission can balance the shortwave absorption.
Without ozone to absorb incoming solar radiation, the temperature does not increase with height.
This simple grey-gas model illustrates this principle very clearly.
Credits¶
The author of this notebook is Brian E. J. Rose, University at Albany.
It was developed in support of ATM 623: Climate Modeling, a graduate-level course in the Department of Atmospheric and Envionmental Sciences, offered in Spring 2015.
%install_ext http://raw.github.com/jrjohansson/version_information/master/version_information.py
%load_ext version_information
%version_information numpy, climlab